Counterion Condensation in Strong, Flexible Polyelectrolytes 



o 

Q 

(N 



O 



I 

o 
o 



> 



I 

O 

o 



X 



Mark J. Stevens and Steve Plimpton 
P.O. Box 5800, MS 1111, Sandia National Laboratories, Albuquerque, NM 87185 

(November 19, 1996) 

We present results of molecular dynamics simulations of strong, flexible polyelectrolyte chains in 
solution with added salt. The effect of added salt on the polyelectrolyte chain structure is fully 
treated for the first time as a function of polymer density. Systems above and below the Manning 
condensation limit are studied. The chain contraction due to added salt is weaker than expected 
from simple screening arguments. The chain structure is intimately tied to the ion density near the 
chain even for chains in the counterion condensation regime. In contradiction to Manning counterion 
condensation theory, the ion density near the polymer chain depends on the amount of added salt, 
and above the condensation limit the chains significantly contract due to added salt. 



Polyelectrolytes are a very important class of poly- 
mers, because they are one set of water-soluble polymers. 
They comprise a class of biopolymers, e.g. DNA, RNA, 
and polysaccharides. Their water solubility makes poly- 
electrolytes an important class of synthetic, commercial 
polymers used for example, as the key water absorbing 
ingredient in disposable diapers. Understanding the ef- 
fect of added salt on polyelectrolyte systems is a basic 
issue for polyelectrolytes. The addition of salt is a key 
means to alter polyelectrolyte structure and properties, 
and distinguishes polyelectrolytes from neutral polymers. 
Nature takes advantage of this in biological systems, and 
the influence of salt is common in industrial applications. 
However, our understanding of polyelectrolyte structure 
is limited, and consequently, so is our understanding of 
polyelectrolyte system properties [^J|,^,||. Recently, 
simulations have been able to determine the structure of 
strong, flexible polyelectrolytes in salt-free solution 
In this Letter, we present results of polyelectrolyte sim- 
ulations with added salt. In particular we examine the 
influence of the ionic distribution on the chain structure. 

Our understanding of polyelectrolytes is poor because 
of the difflculties these systems present to both experi- 
ment and theory. Direct measurements of chain struc- 
ture, particularly at dilute concentrations have yet to be 
done. Recent simulations of salt-free polyelectrolyte sys- 
tems have overcome some of the major theoretical diffi- 
culties 1^. The two main theories Q of polyelectrolyte 
structure have been shown by recent simulations [^|J^] 
to be incorrect for flexible polyelectrolytes. Theoretical 
works have tended to neglect entropy, in part because 
they focused on DNA which is intrinsically very stiff. For 
stiff chains entropy is a small contribution to the free en- 
ergy. In contrast, for flexible polyelectrolytes, treating 
entropy along with the Coulomb interactions is essen- 
tial and has only been done properly in simulation [^,^ , 
although recent self-consistent field theory calculations 
that treat entropy are very promising . 

One of the major theoretical difficulties is the calcu- 
lation of the ionic density about the chain. All calcu- 
lations of chain structure use the Debye-Hiickel (DH) 
approximation. This approximation is a linearization 



of the Poisson-Boltzmann (PB) equation, valid when 
the Coulomb interaction energy is much less than ksT. 
Yet, polyelectrolytes are often highly charged and the 
Coulomb energies can be stronger than fc^T. The PB 
equation [^,p|jic|] is a mean field approximation that ne- 
glects ion correlations, and is valid when typical coun- 
terion separations are larger than the Bjerrum length, 
A = e^/eksT, where e is the dielectric constant of the 
solution (water) and e is the electron charge. More so- 
phisticated methods (e.g. hypernetted-chain (HNC) and 
local density functional theory), have been applied, but 
like the PB method only the problems with simple, fixed 
geometries can be solved pO|JTT| , [T^ . In particular, since 
DNA has a large intrinsic persistence length of a few 
hundred A many calculations treating DNA have been 
performed for cylinders [ pdp^jr^ , p^ . 

Even though it is known that the DH and PB ap- 
proximations can fail, these two approximations pervade 
the field because of their simplicity. Particularly in bio- 
physics, counterion condensation (CC) theory ||l5| and 
the PB theory 0,|l^,|l|l are the two main theories. The 
basis of CC theory is the DH solution of a fixed, lin- 
ear charge distribution (charged line, charged cylinder, 
or a line of point charges) jl^. In Manning's theory, CC 
occurs when ^ = A/a > 1, where a is the charge separa- 
tion distance on the polyelectrolyte. For £, > 1, counte- 
rions condense onto the chain renormalizing the charge 
per unit length so that effective distance between charges 
is o* = A. The condensed ions are contained in a cylin- 
drical volume of radius, Rcc which is salt-independent. 
There is much controversy over the success of CC theory 
12 



19|. Solutions of the PB |]J and HNC equations 
1 11 1 for charged cylinders contradict CC theory showing 
that the ion condensation is dependent on salt concen- 
tration. This is consistent with the CC regime occuring 
where the DH approximation breaks dow n |10|. However, 
CC theory and experiment often agree |18| , |19[ . By sim- 
ulation we can treat the Coulomb interactions without 
approximation and resolve some of these issues. 

Through molecular dynamics simulations, one of us 
has successfully characterized strong, fiexible polyelec- 
trolytes in salt-free solution [ppQ]. The fundamental 
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model of polyelectrolytes in solution consists of chains 
of N charged monomers with N counterions per chain 
and some added salt. The chains have bond length b and 
charge separation a along the backbone. The neutral 
solvent is treated as a continuous, dielectric background. 
The system properties will depend on the polymer den- 
sity, pp, or equivalently, the monomer density, = N pp, 
and the added salt density, ps- Using MD we solve this 
fundamental problem exactly. The results of ealier salt- 
free simulations reproduced experimental measurements 
of the osmotic pressure and the peak position in the inter- 
chain structure factor. Here we treat the case of systems 
with added salt. 

Here, we extend the previous work |^] to include added 
salt H] . The salt ions are modeled the same as the coun- 
terions. In fact, the salt ions oppositely charged to the 
monomers are identical to the counterions. All ions inter- 
act via a purely repulsive Lennard- Jones (RLJ) potential 
with the cutoff at This gives the particles a di- 

ameter of about <T. All ions including the monomers are 
monovalent in this work. The number of added salt ions 
ranges from to 8iVc, where Nc is the number of coun- 
terions which is typically 256 or 512. 

The polymers are modeled as freely-jointed bead- 
spring chains. All beads are charged {b = a). The bond 
potential is the standard FENE (nonlinear spring) po- 
tential with spring constant k = 7e/(T^, and maximum 
extent, i?o = 2cr, where here, as throughout the letter, 
Lennard- Jones units are used The average bond 

length is 6 = l.la. Chains with = 32 and 64 monomers 
are treated. The number of chains is usually 8 with some 
simulations of up to 32 chains. 

The Coulomb coupling strength is determined by the 
Bjerrum length. One set of simulations were performed 
for the same parameter set as in the earlier salt-free work, 
A = 0.83cr In order to study the effects of counte- 
rion condensation, we also performed simulations with 
A = 3.2(7 > a which corresponds to a fully charged chain 
when mapped to sodium polystyrene sulfonate |3|. For 
these simulations we used the particle-particle particle- 
mesh algorithm to calculate the long range Coulomb in- 
teractions mjl^. This algorithm's computation time 
scales roughly as AflogAf, where M is the number of 
charged particles. Such scaling is essential to perform 
these calculations as added salt increases the number 
of charged particles significantly. The simulations per- 
formed here involve up to 16 times the number of charged 
particles as for the corresponding salt-free simulations . 
Even with the almost linear scaling these simulations are 
at the edge of computational capability. 

The temperature is maintained at T = 1.2e, using the 
Langevin thermostat with damping constant T = Ir"^. 
The timestep is 0.015r |2l[. The number of timesteps for 
A^ = 32 and A = 64 is about 3 • 10^ and 8 • 10^ respec- 
tively. This is sufficient to yield reasonable statistics. 

The basic effect of the added salt is to screen the 
charged interactions and contract the polyelectrolyte 
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FIG. 1. The ratio, B, of the end-to-end distance, R, at 
salt density ps to R at no salt for chain lengths A*' = 32 
(open points) and 64 (solid points). Part (a) shows A = 0.83(t 
and (b) has A = 3.2(j. The point types define the monomer 
densities: p^a^ = 0.001(A), 0.01 (□) and 0.02 (Q)- The 
straight line is a guide to the eye and is identical in (a) and 
(b). Error bars are only shown for TV = 32 to maintain clarity. 

chain. We quantified this contraction as a function p^, 
A^ and pm or pp by calculating the average, end-to- 
end distance, R. We calculated the dimensionless ratio, 
B = R{ps)/ R{0), which is the ratio of i? at salt density ps 
to the salt-free value of R. Simulations were performed 
at three monomer densities pm = 0.001, 0.01 and 0.02 
a^'^. For A = 32 all these densities are below the over- 
lap density. For A^ = 64, the salt-free overlap density is 
about Pm = 0.015f7-^ |1. 

Figure ^ shows B plotted against the ratio -y/ Ps/pp- 
Plotting in this manner yields not only a rather linear 
dependence, but also little if any dependence on A^ or 
A within the uncertainty of the data. The square root 
dependence is related to Debye screening, but the ratio 
is not between Debye lengths. The numerator is related 
to the inverse Debye length for the added salt system. 
As — 1/ y'SirXps- However, the denominator is not the 
salt-free Debye length, Ao = l/\/47rAp„i. To obtain a 
quantity proportional to Aq/As, we would need to use 
Pm in Ao instead of pp {pm = Npp). Using p„ yields 
a small, but noticeable A^-dependence. We are skeptical 
that the apparent scaling in Fig. 1 will hold at larger 
A^, and view the data as the best approximate scaling. 
But, as we will see, the screening length is different from 
these simple Debye lengths, and recent calculations show 
the need of a wavelength dependent Debye length . At 
large ps the scaling must break down, because there will 
be a A-dependent crossover to a saturation region. 

A 10% reduction in B occurs at about \J PsI Pp — 8 
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or ps = QApp. For N = 64, this is where ps = Pm] for 
N = 32, Ps = 2pm- To get an appreciable contraction 
beyond the zero salt state, one expects that the screening 
due to the added salt must be stronger than the counte- 
rion screening. In other words, « Aq. Equality occurs 
at ps = Pml'2 which is near where B — 0.9. Thus, when 
the added salt screens the monomer Coulomb repulsion 
about equally to the counterion screening, contraction 
beyond the salt-free state occurs. One would like to un- 
derstand the N dependence better as well as the depen- 
dence on pp that appears important, but simulations at 
larger N are presently too expensive. 

At a sufRciently large ps^ the Coulomb interaction 
should be completely screened and the chain structure 
should be similar to a neutral chain. In terms of De- 
bye screening the monomer-monomer repulsion should be 
completely screened at Ag = a. If the monomer Coulomb 
interactions were treated at the DH level, then at As = a 
the Coulomb interaction becomes just a short-ranged ex- 
cluded volume interaction. In this case, the chains would 
basically be neutral, and a further increase of ps would 
have negligible effect. For A = 0.83(7, this occurs at 
Ps = 0.04cr~'^. At = 32, we have data spanning this 
regime in Fig. 1 at p,„ = 0.02o'^'^ for ps = 0, 2, 4 and 
8pm- No sign of a crossover to saturation appears. At 
Ps = 8pm, R = 10.8(7. The value of R for an equivalent 
neutral chain is 9.7cr. We thus expect the saturation to 
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start for this system at ps ^ 



which is quite a bit 



larger than Debye screening suggests. At pm — O.Olcr^'' 
and ps = 0.04cr~'^, we find R = 12.6a which is far from 
the neutral value suggesting that R will continue to de- 
crease for larger ps . This result shows that screening due 
to explicit salt ions is weaker than simple DH screening. 
In simulations of salt-free solutions, we found that screen- 
ing at dilute densities was stronger when using explicit 
ions than for DH interactions gj]. The results above 
show the opposite occurs for screening due to added salt; 
the interactions between the added salt and the chains is 
screened by the counterions stronger than DH predicts. 

We now turn to results in the CC regime at A = 3.2(7. 
In contradiction to CC theory, Fig. 1(b) shows that the 
chains contract similarly to the A = 0.83(7 case. In CC 
theory adding more salt does not effect the environment 
near the chain. This suggests added salt should only have 
a weak effect on chain structure. On the other hand, PB 
calculations have shown ionic density near the macroion 
is dependent on salt density ||ll| , [l^ . Thus, our results 
agree better with PB predictions. 

We directly examine the ionic distributions in Fig. 
^. The figure shows the normalized radial ion num- 
ber density, n^-^-{r), of all positive ions (i.e. counterion 
and positive salt) about a negatively charged monomer 
for N — 32, A = 3.2a and pm = O.OOlcr""^ at several 
salt densities. This distribution is calculated from the 
monomer-counterion radial distribution function (rdf). 
Since the counterions and positive salt ions are identical 
particles, their rdf's are identical (this was confirmed). 




FIG. 2. Ratio of the positive ion number density about a 
monomer to the bulk counterion number density for A = 32 
at p„ = 0.001(7-^ and A = 3.2a for ps/pm = (Q), 1 (A), 
and 8 (□). This system is in the counterion condensation 
limit and has a dilute polymer concentration. 



The number density is calculated by multiplying the rdf 
by the number of positive ions (the chains are negatively 
charged) . We then normalized by the counterion density. 

The peak heights at different ps values are not iden- 
tical. For Ps — Pc, the peak is about 10% higher than 
the salt-free case. This is not a large difference in peak 
heights. For larger ps the relative peak height increase 
is smaller. However, the peak position is not the only 
relevant spatial regime. As r increases to distances cor- 
responding to just a few monomers away the distributions 
deviate as they must tend toward their asymptotic limit 
of 1 -t- Ps/ Pc- The limit is 9 for the largest ps in Fig. 2. 
At r = 3(7, there is a 50% increase in the ion density for 
Ps — 8pm over ps = 0. Thus, within the volume of the 
chain the ion density increases as Ps increases. Even for 
^ > 1, the repulsion between monomers distant along the 
chain is increasingly screened by increasing ps, and thus 
the chains contract. This contradicts CC theory which 
predicts no salt dependence of the ion density near the 
chain. 

The number of ions near a chain can be directly 
counted. We define the distance between an ion and the 
chain to be the minimum of the distances between the 
ion and each monomer in the chain. Using this metric, 
we can calculate the number of ions within a distance x 
to X + dx of the chain. Normalizing by N, we obtain the 
fractional number, f{x). Integrating f{x) gives the total 
fractional number of ions within a distance, x, which we 
define as F{x). If each chain on average has N counte- 
rions within xn oi it, then F(xn) = 1- In the case of 
added salt, since the counterions and positive salt ions 
are indistinguishable, we calculate F where an ion is ei- 
ther a counterion or a positive salt ion. When F{x) > 1 
for Ps > 0, this means that the number of counterions 
and positive salt ions within a; of a chain is larger than 
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X (a) 

FIG. 3. Fractional number of the positive ions within x of 
a chain for the same systems as Fig. ^. Ps/pm = (cir- 
cle), 1 (triangle), 2 (hexagon) and 8 (square). The Manning 
condensation fraction for these systems is 0.65. 



the number monomers on the chain. This wiU clearly 
happen for large ps- Figure I shows F( x) for the same 
systems as Fig. ||. Because of steric repulsion, there are 
very few ions within la of the chains. Within 2a there 
is a substantial number of ions. At ps = 0, F{2a) = 0.60 
which is less than the Manning fraction for linear chains, 
1/^ = 0.65. With increasing ps, F{2a) increases to about 
0.82 at ps — 8pm, well above the Manning fraction. Fur- 
thermore, within less than 3a, F reaches 1! Thus, the 
situation is more complicated than CC theory predicts. 

All these results imply that extending counterion con- 
densation theory to strong, flexible polyelectrolytes fails 
to predict the correct salt dependence of the chain struc- 
ture. The weakness of the Debye-Hiickel approximation 
is especially strong in the condensation regime and prob- 
ably causes the CC predictions to deviate numerically 
from "exact" values determined by simulation, but it also 
is essential to treat the flexibility of the chain. Given 
our strong disagreement with CC theory, why is CC the- 
ory thought to be well confirmed by experiment? First, 
counterion condensation in some sense clearly occurs near 
^ = 1 in the simulations . It is the statistical physics of 
the counterion condensation which we find to be differ- 
ent for predictions of Manning's CC theory. Part of this 
disgreement occurs because for flexible chains counterion 
condensation alters the chain structure, something that 
beyond the original formulation of CC theory. Measure- 
ment of the polyelectrolyte structure at dilute concentra- 
tions has yet to be done, let alone a measurement of the 
ion distributions near the chain where we find the major 
discrepancies. Much of the agreement with experiment is 
based on the presence of counterion condensation, not on 
the details of CC theory. Furthermore, many experimen- 
tal measurables depend only on the ion distribution in 
the interstitial regime between macroions where the DH 
approximation is often accurate [ p^ , p^ . These quantities 



are often not related to the chain structure or the ionic 
distributions near the chain, and thus are not good tests 
of CC theory's fundamentals. 
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